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Neutral models for the dynamics of a system of competing species are used, nowadays, to describe 
a wide variety of empirical communities. These models are used in many situations, ranging from 
population genetics and ecological biodiversity to macroevolution and cancer tumors. One of the 
main issues discussed within this framework is the relationships between the abundance of a species 
and its age. Here we provide a comprehensive analysis of the age-abundance relationships for fixed- 
size and growing communities. Explicit formulas for the average and the most likely age of a species 
with abundance n are given, together with the full probability distribution function. We further 
discuss the universality of these results and their applicability to the tropical forest community. 
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INTRODUCTION 

An understanding of the mechanisms underlying the composition of a community of competing species is one of the 
important problems of contemporary science. Ecologists would like to understand what determines the stability of 
high-diversity assemblages like freshwater plankton mm or trees in the tropical forest [5] , geneticists are interested in 
the forces that govern genetic diversity within a population, other researchers would like to explain the dynamics of 
different clones within a malignant tumor [4] or the patterns of macroevolution . All these systems are characterized 
by competition between many ’’species” (clones, haplotype etc.) with an underlying birth-death process that allows for 
abundance variations and is combined with the possibility of mutation/speciation by which new species are originated. 

Traditionally, the attempts to explain the dynamics of such communities were based on ’’niche” approach, i.e., on 
the identification of a special feature of any species (like better efficiency in the consumption of a certain resource 
etc.), a feature that provides an advantage for this species and limits the strength of its competition with other 
species [6j. This approach poses many conceptual problems, like the identification of so many different niches mm 
and the stability of a complex community with many competing species [5]. These problems, together with the 
practical difficulty in the calibration of parameters for high-dimensional models, [5], have motivated the search for an 
alternative paradigm. 

During the last decades a radical alternative to the concept of niche-based analysis was presented in various fields. 
Neutral theories [mmi] suggest a new model, under which (at least as a first order approximation) all individuals in a 
community are competing equally for a single resource (e.g., space). To be more specific, there is no niche separation 
(the competition between any two individuals is the same, independent of their species identity) and the fitness of 
all individuals is, again, the same (demographic equivalence). In a neutral game, the outcome, like the abundance of 
a species or its chance to perish, depends solely on the stochastic processes in the system, and in particular on the 
demographic fluctuations that characterize the birth-death-mutation dynamics. 

An important distinction between various types of neutral theories has to do with the overall size of the community. 
In ecology and genetics the research is focused on communities that have, more or less, a fixed size, assuming that the 
overall biomass of the community have already reached a stable equilibrium. On the other hand, when the neutral 
theory is applied to macroevolution [^, to the dynamics of networks or to surnames in human population 

[141 - IH)j . one has in mind an exponentially growing community. Recently, a few authors applied the neutral theory 
to a growing malignant tumor, trying to understand intra-tumor heterogeneity and mutant allele frequencies [4]. In 
particular, the heterogenous response of individuals to drug therapy as was interpreted as reflecting the abundance 
of resistent cells, assuming that sensitive and resistant clones have the same demographic parameters before the 
treatment [IIHIS]- 

In many cases the success of the neutral theory in explaining community composition (especially species abundance 
distribution (SAD)) is quite impressive, given its minimalistic assumptions and the small number of parameters 
involved [151120j . However, the predictions of the neutral theory should be examined also in relation to the dynamics 
of empirical systems m- In this context, one of the main differences between neutral and niche theories has to do 
with the relationships between the age of a species (the number of generations elapsed since the speciation event, 
i.e., the birth of the founder) and its abundance. In niche theories there is a very weak correlation between these 
two quantities: the model supports an attractive fixed point that depends on the relative fitness of the competing 
species, and the equilibration time is almost independent of the number of individuals. Under neutral dynamics, on 
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the other hand, the system is governed by stochasticity. Accordingly, in a fixed size community the time needed for 
a single founder to establish population of size N is typically N generations, and in a growing community the age 
of a species is positively correlated with its age. Indeed, one of the known difficulties of Hubbell’s neutral theory 
of biodiversity, has to do with the age-abundance relationships, since the N generations requirement appears to be 
unrealistically long [52H23]. For growing communities, like in the case of macroevolution, one can find contradicting 
statements about the correlation between taxon’s age and its size [3 [13 [in]- 

This paper is dedicated to the age-abundance problem under neutral dynamics. We extend the results of a recent 
work by Chisholm and O’Dwyer |24j . who provide an analytic expressions for these relationships in a fixed size 
community. Here we consider also growing (or shrinking) communities, showing that under neutral dynamics, the 
species abundance distribution determines unequivocally the relationships between abundance and mean age of a 
species. We further explore the conditions under which these results are universal, i.e., independent of the details of 
the process (like overlapping/non-overlapping generations, or the statistics of offspring per individual). 

This paper is organized as follows: in the first few sections we present a solvable model for Wright-Fisher dynamics 
(non-overlapping generations) with geometric distribution of offspring. , the chance of a species to have abundance 
n, s generations after its origination is calculated, along with the maximum likelihood expression for the age of the 
species, s, given n. We then proceed to calculate T*{n), the average age of a species given its abundance, showing 
that this quantity is determined by the species abundance distribution. Finally, we discuss the universality of our 
results and its applicability to the empirical findings in the Amazon forest. 


A SOLVABLE MODEL 

The master equation for a continuous time neutral dynamics of a single population with pure demographic noise 
(a Moran process) is given by, 


=-2mPjn +{m + l)Pm+i + (1) 
at 

where Pm is the chance that the species is represented by m individuals. This equation corresponds to a process in 
which every individual tosses a coin at random, dying if the outcome is tails and producing an offspring when the 
outcome is a heads. The probability (per individual per unit time) for such an event sets the timescale in Eq. Q. 
Accordingly, the chance that an individual produces n offspring during its lifetime (from birth to death) is given by 
a geometric distribution, Pn'^ = 1/2”+^, where the superscript (1) indicates that this is the probability to produce n 
offspring during a single ” generation”. 

Here we consider a variant of this model, with the same offspring statistics but with non-overlapping generations 
(Wright-Fisher dynamics). In this version at every generation each individual produces n offspring with a chance 
Pn = 2“^"+^^ and dies. The only difference between the continuous and the discrete time versions is the definition of 
a generation, and as we shall see below in the relevant regime of parameters this modification has no importance. 

For the geometric distribution it is quite easy to calculate the one generation generating function. 


OO ^ 

GW(a:) = ^a:'^pW = ^-, (2) 

2 — X 

n—O 

and the structure of this generating function [27j allows for an explicit formula for the generating function after s 
generations, which is obtained by iterating s times. 


G^^\x) = =G(G(G...s times (G(x)))) (3) 

71 — 0 

s — {s — l)a: 

(s -I- 1) — sx 

(s) 

The probability of a single individual to produce n offspring after s generations, Pn , is the coefficient of x" in the 
series expansion of G^^^x). This yields (for n > 0), 
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while for n = 0, = s/(s + 1). Accordingly, the probability of survival s generations decays in long times like 1/s, 

in agreement with the general results of Galton-Watson theory [28] . 

To generalize this model we consider a case where the average number of offspring produced by a single individual 
is 1 + 7 (so 7 is the growth/decay rate of the population); the chance of success in the corresponding Bernoulli trial 
should be z = (7 + l )/(7 + 2). In addition we consider a process with mutation/speciation events: the chance that 
an offspring does not belong to the same taxon of its mother is v (no recurrent mutations, and mutant is the founder 
of a new species). For this birth-death-mutation process the average number of non-mutant offspring per individual 
is Po = 1 + r = (1 -I- 7)(1 - ly). 

For the Bernoulli process considered here the probability for n non-mutant offspring is given by: 


p(i) 

n 




(l + Po)-+l’ 


(5) 


hence the generating function is. 


G^^\x) 


1 

1 -I- Po — RoX 


( 6 ) 


Although the recursion relation is less trivial, one can still derive the generating function after s generations. To do 
that, it is useful to notice that during the recursion process (plugging instead of a; in (|^ to get G^^^ and so on) 
G retains its general form, 


^ _as+ PsX 
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and the recursion from the s to the s -I- 1 generation satisfies. 
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with the initial conditions oi = 1, /3i = 0, 71 = 1 -|- Ro and i5i = —Ro- Accordingly, 


(7) 


( 8 ) 


^ (Pg - 1) - (Pg - Ro)x 
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and. 
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( 10 ) 


Again, the last result is valid for n > 1, while for n = 0 one obtains Pg*^ = [1 ~ 7?g]/[l — Pg"*"^]. 


ABUNDANCE OF A CLONE AS A FUNCTION OF ITS AGE 

Armed with the formula obtained in the last section, one can get some basic intuition regarding the age-abundance 
relationships. The average abundance of a species s generation after point speciation (conditioned on non-extinction, 
i.e., the average number of descendants of a single individual after s generation, given that this individual has at least 
one living descendant at this time) is, 


n(s) 
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r<0 r=0 r>0 



FIG. 1: The average size of a (non-extinct) species s generation after its origination, n(s), as a function of s. ... present the 
results of a numerical simulation, while the full line is the analytic expression, Eq. [m In the left panel, f?o < 1 and r < 0 
meaning that, although the overall size of the community is growing 7 = 10 “^, the abundance of every species is shrinking 
on average as the demographic losses due to mutations 1 / = 2 * 10~^ are stronger. In that case, even species that survive 
for a long time cannot cross (on average) the abundance l/|r|. In the middle panel the behavior of h(s) marginal case r = 0 
(7 = 10“^, u = 0.99*10“^) is graphed, where now the growth is strictly linear, while in the left panel (7 = 2*10“^, v = 0.0102) 
the crossover to exponential growth happens, again, at 1/r. Overall, the growth is controlled by the ’’deterministic” parameter 
r after l/|r| generations, and before this stage the growth is linear, independent of r, and is related to the purely stochastic 
dynamics at r = 0. The numerical results presented here were obtained by simulation a Wright-Fisher process with the 
appropriate values of 7 and v, starting from a single individual and averaging over 10 ® trials. 


To understand the implications of (111, let us start with the marginal case i?o = 1 (or r = 0). In this case the 
overall community is growing (7 > 0) but the average size of a single clone is kept fixed on average, as the effect of 
growth is balanced by mutations. In this limit the average abundance of a surviving clone is growing linearly with 
its age s, n{s) = s + 1. On the other hand, the average abundance of all clones (including those who went extinct 
and has abundance zero) is time independent, since the chance of surviving until s scales like 1/s. This feature is 
demonstrated in the middle panel of Fig[2 


The marginal case Ro = 1 is, of course, special and requires a fine-tuning balance between mutations and growth 
rates. Still, in many common situations the relevant values of both 7 and v are much smaller than unity (see, in 
example, the cases considered in m). hence r << 1. In these cases, as long as |r|s << 1 the surviving clones are 
still growing linearly in time as in the marginal case. This linear pattern breaks down for |r|s>l: ifr>0 the clone 
switches to exponential growth, n{s) ~ exp(rs) (Fig. right), while for r < 0, n{s) saturates to 1/r at large s, as 
demonstrated in the left panel of Fig. This latter case includes the scenario considered in Hubbell-Kimura neutral 
theory, i.e., a fixed size community, 7 = 0 , and finite chance of mutation, ly > 0. 


AGE-ABUNDANCE RELATIONSHIPS: A MAXIMUM LIKELIHOOD APPROACH 


As seen in the last section, our solvable model allows us to calculate the average abundance of a species given its 
age. However, many practical problems have to do with the opposite question: given the abundance of a species or a 
clone and its demographic features represented by r, what is the best estimate one may suggest for its age. 


Perhaps the simplest approach towards this problem is to implement Eq. (101 and to calculate, given r, the value 
of sml that gives the maximum likelihood for n. It turns out that. 


sml — 


n/2 

ln{n\r\) 


hi 


n|r| < 1 
nlrl > 1 


( 12 ) 


Note that the expression for Sml is independent of the sign of r, i.e., it takes the same value for growing and shrinking 
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clones. This feature is obvious in the linear regime n\r\ <C 1, where the system still behaves as if it is at criticality. 
Quite surprisingly, however, even when n|r| ^ 1 sml is still independent of the sign of r despite the pronounced 
difference in the actual behavior of the two systems. When r > 0 the long term growth is exponential, so the time 
to reach abundance n is logarithmic in n. When r is negative, as demonstrated in Fig clones with abundance 
n > l/|r| are exponentially rare and do not contribute to n. Still, the most likely course of events that leads to this 
rare scenario is a rapid exponential growth, so the age-abundance relationships are still logarithmic. 

For the special case of a fixed size community (7 = 0 , i/ > 0, r = —v < 0) the species abundance distribution 
(SAD) is known to be the Fisher log-series exp{—vn)/n. Accordingly, species with n > 1/v are exponentially rare, and 
typically the range of possible sizes of species is limited by l/|r|. As followed from ( |l^ , if such a rare event occurs 
the age of the rare species scales with the logarithm of its abundance (as opposed to the linear scaling for species in 
the bulk of the SAD). The relatively young age of the exceptionally abundant species cannot solve the species-age 
problem, since the chance to observe species with n >> 1 /iz is tiny. 


THE AVERAGE AGE OF SPECIES OF A CERTAIN ABUNDANCE 


In this section we would like to calculate the average age of a species of abundance n (as opposed to the maximum 
likelihood method implemented in the last section). When the distribution P{s\n) is sharply peaked, the mean and 
the maximum likelihood are almost identical. On the other hand, when the distribution is wide and skewed the 
average provides additional information and facilitates the estimation based on finite number of data points. 

To begin, we assume that the number of mutations/speciations per generation is given by i/J(s), where J is the 
overall size of the community. When the size of the community is growing exponentially one may assume further that 
J(s) = exp(js). Given that, the expression P{n\s) (the chance that a species, observed at a size s, was originated n 
generations ago) may be obtained from Pn'^ using elementary Bayesian inference |24) . 


P{s\n) 




2^s=i 


(13) 


The denominator in this expression is proportional to the average number of clones/species of size n at the s generation, 
a quantity that was already calculated in the context of species abundance distribution for growing and shrinking 
populations [HIlsllISllIl]. Since Pn'^ is known, P{s\n) may be written explicitly, and we will provide the appropriate 
expression at the end of this section. For the moment let us calculate analytically the average lifetime of a species 
given its abundance as suggested in [ 21 ], 


T*{n) = sP(s|n). 


(14) 


Implementing a fictitious parameter 7 , T* (n) may be written in a simpler form (this step is analogous to the 
derivation of the free energy from the partition function in equilibrium statistical mechanics), 


T*{n) = 




sPk^^e- 




pk"'‘e-y^ 



c-1 


E 




(15) 


Note that Pkf'^ is also 7 dependent, as shown in the last section, via r = (1 -|- 7)(1 — v) (or, if both rates are much 
smaller than one, r « 7 — i^). To get the expression in the r.h.s. it thus crucial to make a technical distinction between 
the growth rate of the whole community 7 and the growth rate of a single clone, j — 1 / which is hidden in h As 
will be soon explained, this technical point facilitates the analytic calculation. 

As mentioned, the argument of the log in the r.h.s. of (15) is, up to constants, M„, the average number of species 
with abundance n, which is simply the number of new species (mutants) at a certain generation [u times the number 
of birth which is the size of the community at this generation) multiplied by the chance that a mutant will generate 
a family of size n when the community is sampled. 


t-i 


= y ve 


S =1 


7(t-s)p(s)_ 


As shown in [HdllllS], this quantity satisfies. 
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When the SAD reaches a stable distribution, the chance of a clone, picked at random (independent of its abundance) 
from the list of living clones, to have abundance n is time-independent (this, of course, cannot be true for the wild 
type as will be discussed soon). Accordingly, for a growing community the number of n-clones, M„, must grow 
exponentially like expifjt). Hence, 


yielding dansnii]. 


2 ^ 
^M{n) = (nM(n)) + (^ “ 7) 


n \ It — cr 


(18) 


(19) 


where U is the Kummer function [^ . 

Before plugging this expression into (15), 
condition is, 


we should consider the normalization factor A. 


The normalization 


nM{n,s)dn= / nM{n)e^'^dn = J„{s), 


( 20 ) 


where Ju{s) is the overall abundance of all the mutant species, i.e., the size of the community apart from the wild 
type, 


Y(s) = [ = 
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7S 
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Combining (20) and (21) one obtains: 


( 21 ) 


nM{n)e^^ dn = 


Plugging M(n) into this expression one can find A, yielding. 


M{n) = ( 1 + 


|7- 
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7 - 7 -I- 
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Finally, the average age of species with abundance n is given by: 
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l7- 
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\l- 




( 22 ) 


(23) 


(24) 


To compare with the maximum likelihood results presented above one would like to translate 'y — v to r. Note, 
however, that the Fokker-Planck expression holds only for small values of 7 and u, i.e., when is negligibly small. 

As above, we will consider two limits, |r|n << 1 and |r|n >> 1 . Although the results are relevant for both 
positive and negative r, our interest is mainly in the positive r case so we will explain the intuitive argument in the 
corresponding language. In the first limit rn << 1 we deal with clones that are still in their linear growth phase. 
Expanding T* for small rn one obtains: 


T*{n) « —n I log(nr) + — + f —'j -|- 


7 


■jijjA) ^ 


Htt'i 


+ 27_e - 1 , 


(25) 


so the leading term scales like nln(n), as opposed to the linear dependence obtained by maximum likelihood. 
It is interesting to note that, in the limit of fixed size community (7 = 0, r = v), 


(26) 


T*{n) « -n(log(rn^) + - I), 

and for large n the identity « 7 ^; -|- ln{n) leads to the result of [23]: 

T*{n) = n [I — Hn — In {v)]. 


(27) 
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In the large rn limit the species is in its exponential growth phase so the age should be logarithmic in the abundance. 
Indeed, for large rn the Kummer function may be approximated by its power law asymptote to yield, 


log(nr) 

r 


—H:t + 7 e 
r 


(28) 


and the leading term is identical with the maximum likelihood estimation. Since our theory depends only on the 
absolute value of r one may take the limit 7 = 0 , \r\ = n to obtain, again, the result of 


T*{n) K. - (log(i/n) + 7 b) . 


(29) 


Finally, since we have an expression for Pn \ the the probability function P{s\n) may be derived from (13). For 
example, the geometric model yields. 


P{s\n) 


T-2n I " 
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(30) 


DISCUSSION 


The main goal of this work was to extend the work of [23] to communities with finite growth (or decay) rates. 
To do that we used two approaches: the first was based on a solvable model, the Wright-Fisher dynamics with 
geometric distribution of offspring per individual. This model allow s fo r an exact solution for the probability of a 
single individual to produce n offspring after s generations, Pn'^ (Eq. 10). The emerging picture is that the surviving 
species grow linearly in time as long as s < l/|r|; at longer times they either grow exponentially (if r > 0 ) or saturate 
(if r < 0). Eq. (10) provides also a maximum likelihood estimation of the age of a species given its abundance. 
Quite surprisingly it turns out that the result does not depend on the sign of r. This apparently reflects some hidden 
symmetry in this system. 

Our second approach was based on the Bayesian method suggested in |24|, extending it to the case of growing and 
shrinking populations. In this case we have calculated T*(n), the average age of a species with abundance n. To do 
that we have implemented a trick, based on the distinction between the growth rate of a population 7 and the growth 
rate of the community 7 . With this trick, the problem looks very much like the calculation of expectation values for 
a statistical physics system at equilibrium, and the solution depends only on the species abundance distribution of 
the model. For growing or shrinking communities the SAD has already been calculated [131 113 IE] : and this led us 


to the expression (24) and its low/high abundance asymptotics, including the results of [24] as a special case. 
At this point, we would like to add two technical comments. 


1. First, while the first approach was based on a specific model (geometric distribution for the number of offspring) 
for which the average number of same species offspring, i?o = (1 + 7)(1 —i^), determines the variance Ro{1 + Rq), 
the analysis of the second approach is based on the SAD, given in general terms without restriction to a certain 
relationships between the mean and the variance. The justification for that has to do with the universality of 
the problem, as analyzed in detail in m- Strictly speaking, the use of a Fokker-Planck equation is justified 
only in the limit of small r and Oil) variance such that may be considered as a continuous function of n 
and s. Therefore, given any distribution of offspring P„ with mean 1 + 7 and variance cr^ , on e can implement 
the Fokker-Planck equation to obtain Pn'’ = aP^'^l (an) where (n) fits our result (Eq. 10) and 


Po(l + Rq) 

(1 — n)^cr2 _|_ jy/jp 


This feature is demonstrated in Fig. [^for geometric and Poisson offspring distributions. 

2. Another comment has to do with the applicability of the species abundance distribution for a growing community. 
The SAD that we have analyzed is based on a steady state solution of a Fokker-Planck equation. However, while 
for a fixed size community, when the maximum abundance of a species scales with l/z/, the equilibration time is 
also finite, for a growing community with r > 0 there is no finite limit to the abundance and the equilibration 
time is always too small for the first few species (like the wild type and the first mutants). Accordingly, the 

















FIG. 2: n(s) (conditioned on survival) for geometric (red) and Poisson (blue) distributions of the number of offspring per 
individual. In the inset, the number of offspring as a function of time, averaged over 10* histories, is plotted for two communities. 
Both communities have 7 = 2*10“^ and u = 10“^, and the only difference is the distribution function. Apparently, the geometric 
species grow much faster; the reason is that the chance of extinction for a geometric species is higher, and the surviving species 
abundance should compensate it in order to preserve the fixed mean. In the main panel both histories are plotted again, now 
the Poisson (red) points are divided by a factor of where for Poisson cr^ = (1 — + uRq and = Ro{Ro + 1) is 

the geometric distribution variance. Both datasets collapse, and fit the analytic expression [m 


Kummer-like SAD analyzed here is inappropriate for the most abundant species, like those that were considered 
by [niEO]. Note, however, that the width of the distribution for abundance values for these species is very 
wide, making a reliable inference of the age of the species from its size quite difficult. 

Finally, we would like to discuss the applicability of a model with growth/decay of the population to the species-age 
problem that arises in the context of Hubbell’s neutral model. As we have already seen, the typical time needed for 
a species to reach abundance n is 0{n) generations as long as nr < 1 , which is always the case when 7 = 0 and 
the largest abundance in the system is OiXjv), the regime considered in the neutral theory of biodiversity. If, for 
example, one considers a set of 10® conspecific trees for a frequent species in the Amazon basin (this is close to the 
contemporary figure, see recent survey in 0), with about a 50y generation time, the neutral theory suggests 50 • 10® 
generations, more than the age of the universe. The tempo of the dynamics may be accelerated if one assumes a 
very large value of cr® (as suggested, essentially, in [31], see |24|. Note that the n scaling is renormalized by cr® as we 
discussed in the context of universality) or by keeping the generation time as a free parameter (see, e.g., |32|), but 
any of these approaches carries its own difficulties when compared with realistic timescales and variances in empirical 
systems. 

One of the original aims of this project was to check if one can solve this timescale problem by assuming that the 
trees in the Amazon basin, say, are a growing community. In many models evolution takes place as species increase 
their carrying capacity (K-selection), so at least in principle it is possible that the population was much smaller in 
previous times. However, when we tried to apply our analysis to the estimations suggested by |3] we failed to find an 
appropriate regime of parameters. Specifically, a fit of Kummer function to the observed SAD result in very small 
values of 7 , values that cannot solve the age-abundance problem. 
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